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ABSTRACT 

An  implicit  approximate  factorization  (AF)  algorithm  is  constructed  which  has  the  following 
characteristics. 

•  In  2>D:  The  scheme  is  unconditionally  stable,  has  a  3  x  3  stencil  and  at  steady  state  has  a 
fourth  order  spatial  accuracy.  The  temporal  evolution  is  time  accurate  either  to  1st  or  2nd 
order  through  choice  of  parameter. 

•  In  3-D:  The  scheme  has  almost  the  same  properties  as  in  2-D  except  that  it  is  now  only 
conditionally  stable,  with  the  stability  condition  (the  CFL  number)  being  dependent  on  the 
*cell  aspect  ratios,”  Ay/Az  and  Az/Az.  The  stencil  is  still  compact  and  fourth  order  accuracy 
at  steady  state  is  maintained. 

Numerical  experiments  on  a  2-D  shock-reflection  problem  show  the  expected  improvement  over 
lower  order  schemes,  not  only  in  accuracy  (measured  by  the  Lj  error)  but  also  in  the  dispersion. 

It  is  also  shown  how  the  same  technique  is  immediately  extendable  to  Runge-Kutta  type  schemes 
resulting  in  improved  stabUity  in  addition  to  the  enhanced  accuracy.  ^ 
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1  INTRODUCTION 

It  can  be  shown  [1]  that  numerical  approximations  to  the  linearized  Euler  equations  of  gas  dynamics 
give  rise  to  dispersive  errors  which  in  the  2-D  supersonic  case  depend  on  a  similarity  parameter 
K  =  (Ay/Ax)\/M2o  ~  1  (under  the  assumption  v  C  u  everywhere,  where  u  and  v  are  the  x 
and  y  components  of  the  velocity  vector).  The  difference  between  the  dispersion  relations  of  any 
numerical  algorithm  and  that  of  the  original  partial  differential  equations  can  be  plotted  as  curves 
in  the  Fourier  plane  with  k  as  a  parameter. 

In  particular  the  results  in  [1]  indicate  that  for  central-difference  schemes,  the  dispersive  error 
are  contributed  mostly  by  the  third  power  of  the  errors  in  the  Fourier  variables  $  and  It  is, 
therefore,  natural  to  think  of  fourth  order  spatially  accurate  algorithms  as  having  better  dispersive 
properties.  By  utilizing  the  structure  of  the  Euler  equations  one  can  obtain,  on  a  cartesian  grid,  a 
fourth  order  approximation  which  instead  of  using  a  5  x  5  stencil  (and  5  x  5  x  5  in  3-D)  relies  on 
a  compact  support  of  3  x  3  (and  3  x  3  x  3  in  3-D).  The  advantages  of  the  combination  of  fourth 
order  accuracy  together  with  compact  support  are  quite  obvious  in  terms  of  total  computer  work 
and  memory. 

In  Section  2,  we  describe  the  construction  of  an  approximate  factorization  (AF)  central  differ¬ 
ence  scheme  and  examine  its  theoretical  (linear)  stability  properties.  In  Section  3,  we  derive  the 
corresponding  4-step  Runge-Kutta  scheme  and  show  how  the  Jameson-Schmidt-Turkel  algorithms 
[2]  may  be  easily  modified  to  that  form  which  has,  in  addition  to  the  higher  accuracy,  markedly  en¬ 
hanced  stability  limits.  In  Section  4,  we  describe  some  numerical  experiments  using  the  AF-version. 
Section  5  summarizes  our  findings. 

2  DERIVATION  OF  THE  APPROXIMATE-FACTORIZATION 
SCHEME 

2.1  The  Two-Dimensional  Case 

Consider  a  general  hyperbolic  conservation  law  in  2-D: 


«!  +  /*  +  9t=  0. 


In  the  case  of  Euler  equations,  for  example,  the  vectors  u ,  / 
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where  p,  u,  v,  E,  and  p  are  respectively  the  density,  velocity  components  in  the  x  and  y  directions, 
the  total  energy  per  unit  volume  and  the  pressure.  In  addition  there  is  the  equation  of  state  relating 
algebraically  (in  the  case  of  gas-dynamics)  the  internal  energy  to  the  pressure  and  density.  One 
may  also  write  (1)  in  non-conservation  form  as 

ut  -l-A  u,  +B  Uy=  0 


where  A  and  B  are  the  Jacobians  of  /  and  9  with  respect  to  u  . 


(3) 


Consider  first  a  forward-time  second  order  central  finite  difference  approximation  to  (1) 

At  ^  A*  Ay  ^  ^ 

where  we  have  dropped  the  sup-arrow  indicating  a  vector  and  used  the  conventional  differencing 
notation 

“7*=  «i.k  =  "^0  ■  (5) 

We  have  a  cartesian  grid  with  nodes  at  z  =  yAz,y  =  JbAy  and  t  =  nAt.  We  shall  now  introduce 
the  usual  shift  operations  notation: 

M.U  =  Hy  =  iKt+^  +  (6) 


Eq.  (4)  then  may  be  written  as: 


^v-“y,fc+i 


Dtu  -  -  -^liySyg  (7) 

where  A  =  At/Az  and  Bl  =  Ay/Az.  If  we  take  a  Taylor  series  expansion  about  the  grid  point 
(^iViO  =  ^  construct  the  modified  equation  corresponding  to  (7): 
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Thus  if  we  wut  to  approximate  (1)  to  a  higher  than  second  order  (and  in  particular  to  fourth 
order  in  space-see  comment  in  Introduction  concerning  dispersive  errors)  we  must  modify  (7)  by 
subtracting  out  the  terms  on  the  right  hand  side  of  (8).  At  this  point  we  realize  that  using  a 
straightforward  approach  to  approximating  and  g^yy  will  leiul  to  bigger  stencils.  However, 
using  the  original  partial  differential  equation,  (1),  we  have 

/*«*  =  -«<««  -  gyzz  (9) 

and  similarly 

Since  f^xz  ^nd  gyyy  need  be  approximated  only  to  second  order  (because  of  their  coefficients,  Az^/3 
and  Ay’/3)  the  required  stencil  for  the  terms  on  the  right  hand  side  is  only  3x3.  (This  observation 
was  previously  made,  in  another  context  by  Jones,  et  al.  [3].)  So,  replacing  fxxz  and  gyyy  in  (8) 
by  -S^Dtu/^x*At  -  Slny8yg/^x^^y  and  -8^Dtu/^y^^t  -  6y*/i*5,//Ay*Az,  respectively;  and 
also  replacing  utt  by  D|Ut/At  — *  Dt{-fx  -  gy)/^*  -*  -Dt  [(Mx^b//Az)  +  {t^^yg/^y)]  we  obtain 
the  following  approximation  to  (1)  which  is  spatially  fourth  order  accurate  and  temporally  second 
order  accurate: 
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Using  (3)  and  the  definition  of  Df  we  rearrange  the  above  into  the  delta  form 

[/  +  |(a/..S.  +  ICS.)  +  i  (^5|  +  Bel)]  (“S‘  - 

=  (/  +  yi)  /&  -  ^  (^  +  hn)  »x.  (11) 

where  /  is  the  identity  matrix.  This  is  an  implicit  unfactored  scheme.  Since  we  are  interested  in 
marching  to  steady  state  we  approximate-factor  the  left  hand  side  to  obtain: 

= -A  {»^«.  (/ + i«;)  /,% + iv,  (/ + }  •  (12) 

We  introduced  into  (12)  the  parameter  o'.  If  <r  =  1  then  we  retain  the  temporal  second  order 
accuracy  while  o  =  2  gives  first  order  in  time.  Note  that  (12)  involves  the  inversion  of  block- 
tridiagonal  matrices.  The  right  hand  side  represents  the  steady  state  operator  to  fourth-order.  The 
whole  scheme  involves  only  a  3  stencil. 

2.2  The  3-D  Case 

The  starting  point,  corresponding  to  (1)  is 

««  +  /*  +  +  h»=  0.  (13) 

Following  the  same  steps  as  in  the  2-D  case,  we  get  the  modified  equation 

,  ,  [At  Ax*  Ay*  As* 

“»  +  /*+  5v  +  *»  =  -  +  -^/xx»  +  -^9wv  +  -^hg„  (14) 

where  using  (13),  we  have 

fzxx  —  ~^txx  ~  9yxx  ~  h»xx  (^®) 

9yyy  =  ” Olyy  —  ^xyy  ~  /xyy  (16) 

~  fxXM  ~  9y*f  (17) 

Repeating  all  the  steps  leading  to  (12)  we  obtain  its  three  dimensional  analog: 

(/  +  l-Ml  +  \XA^e.)  (/  +  (/  +  ic«J  +  (.Ji‘  -  .J,.,)  = 
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(18) 

In  (18)  the  matrix  C  is  the  Jacobian  d  /d  u  and  the  shift  operators  are  defined  in  a 
manner  analogous  to  (6).  Q  is  the  cell  aspect  ratio  in  the  t  direction,  Q  =  At/^x. 


3.S  Stability  Analysis  for  the  2-D  Case 

We  consider  the  linearized  (i.e.,  frozen  coeflBcients)  version  of  (12).  We  carry  out  a  Von-Neumann 
axialysis  in  the  usual  manner  by  casting  (12)  in  Fourier  space.  A  typical  Fourier  component  of 
is  given  by 

(19) 

where 
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(20) 

with 

and  — oo  <  i,m  <  oo. 

The  mapping  of  the  various  shift  operations  is  as  follows: 

(21) 
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and  so, 

-*  2ir)y/l  -  17*. 

(25) 

We  also  define  the  amplification  matrix  G  by 

^n+i  _  Qf^n  these  notations. 

the  frozen 

coefficient  version  of  (12)  is  mapped  into  the  Fourier  space  as  follows: 

(/  -  (/  -  |v  +  -  »*)  (G  -  n 

=  -2ix  [Ae,/rii  (i  - 1,»)  +  ,  (26) 

In  the  general  case  the  matrices  A  and  B  do  not  commute,  thus  rendering  the  analysis  of  (26) 
almost  intractable.  It  is  instructive,  however,  to  consider  the  scalar  case.  Since  the  aspect  ratio 
Hi  =  Ay/Az  is  arbitrary,  and  since  the  original  partial  differential  equation  (1),  in  the  scalar  linear 
case,  could  be  transformed  to  the  wave  equation 

U|  +  U|  +  tt  j  =  0 

with  t  =  At,  i  =  X  and  y  =  Ay/ B,  we  can  without  loss  of  generality  (in  this  special  linear  scalar 
case)  rewrite  (26)  as: 

-  {«)  (l  -  (°  -  ‘) 

= -2.A  [{v/r^(i  -  li’) + - !«’)] 

Equation  (27)  can  be  rearranged  to  solve  for  the  amplification  factor  G  : 


(27) 


where 


(29) 


" = ‘  -  |f’  -  k + lev  - 

and  6  is  the  Fourier  map  of  the  steady-state  operator,  i.e.: 

»=2[Aev^r^(i-|,>)+^v^^(i-|e>)].  (30) 

We  see  immediately  from  (28)  that  for  o'  =  1  (i.e.,  the  case  of  second  order  temporal  accuracy)  we 
have  a  Crank-Nicholson  type  scheme  with  ||Cr||  =  1.  For  <r  >  1  (i.e.,  first  order  accuracy  in  time), 
we  have  ||G||  <  1.  We  have  thus  demonstrated  the  unconditional  stability  of  (27)  for  all  values  of 
the  cell  aspect  ratio. 

2.4  Stability  for  the  3-D  Case 
The  three  dimensional  analog  to  (27)  is 

(1  -  +  ioXiyJl  -  f*)(l  -  |ij*  -f-  -  »7*)(1  -  |f’  +  •>^fV^?’)(G  -  1)  = 

=  -KA  [f,/rr?(i  -  fo’  -  |j>) + -  ft*  -  |e’) + ^V^(1  - 1£’  - 1-.-) 

(31 

where  ;  is  the  Fourier  dual  variable  defined  analogously  to  (  and  i;. 

The  stabUity  of  three  dimensional  amplification  factor  G,  as  defined  by  (31),  is  difficult  to 
analyze  and  we  resorted  to  numerical  evaluations  of  jG]*,  using  up  to  8  x  10*  Fourier  modes.  The 
numerical  study  of  (31)  was  carried  out  on  the  Cray  2.  We  found  that  as  in  the  case  of  forward 
Euler  approximate  factorization  second  order  scheme  [4],  the  amplification  factor  was  conditionally 
stable.  For  example  for  IR  =  Q  =  1,  the  stability  limit  is  A  <  .43.  These  stability  properties  are 
obtained  with  the  aid  of  artificial  viscosity  (AV)  term  which  is  of  sixth  order  but  still  resides  on 
the  compact  stencil.  The  AV  term  is  added  to  the  right  hand  side  (explicit  term)  of  (18)  and  is  of 
the  form 

(32) 

where  is  of  order  unity.  We  found  1.5  <  <  2  to  be  most  efficacious.  Without  the  artificial 

viscosity  term,  the  allowed  value  of  A  is  about  one  order  of  magnitude  less  (e.g.,  for  i2  =  Q  =  1,  A 
without  using  AV  is  about  .035). 


3  COMPACT  FOUR  STEP  FOURTH  ORDER  RUNGE-KUTTA 
ALGORITHM 

The  basic  four  step  explicit  Runge-Kutta  scheme,  as  proposed  in  [2],  takes  the  form: 

„(i)  =  „(")  _  1 

=  ut")  — 

(33) 

u(»)  =  ot")  -  ^AtR{u(*)) 
u(<)  =  „(«)  _  AtR(u(*)) 
u("+i)  =  u(^) 


where  R  is  the  finite  difference  (or  finite  volume)  representation  of  the  steady  operator.  For  example, 
in  the  2-D  second  order  spatial  accuracy  case 


(34) 


where  /  and  g  are  the  fiux  vectors  encountered  in  section  2.  If  we  want  fourth  order  spatial  accuracy, 
then  it  follows  directly  that  the  residual  R  is  modified  to  read 


R  = 


Hz^x 


(i+k’)/+ 


Ax  '•  '  6 


(35) 


and  we  still  retain  the  compact  support. 

Similarly,  in  the  three  dimensional  case  we  have 


(36) 


Thus  (33),  with  R  given  either  by  (35)  for  the  2-D  case  or  by  (36)  for  the  3-D  case,  retains  all  the 
features  of  the  second  order  scheme  but  gains  us  the  fourth  order  accuracy.  In  addition  one  can 
easily  verify  by  simple  analysis  that  for  a  given  cartesian  grid  and  flow  conditions  the  new  fourth 
order  formulation  enhances  the  stability  condition.  In  the  2-D  case  we  have,  using  (35)  rather  than 
(S4) 

it-!"*”*" pn 

(^0(2)  (^*)2nd  order 

In  the  3-D  case  the  gain  is  even  more  favorable, 


(At)4 

(At), 


=  1.66. 


(38) 


Thus  the  algorithm  efficiency  gains  are  two  fold.  First,  for  a  given  acceptable  error  level  the  fourth 
order  accuracy  allows  a  coarser  grid,  i.e.,  fewer  node  points.  Second,  not  only  At  is  increased  due 
to  the  larger  cell  size  but  in  addition  it  gains  due  to  (37)  (or  (38)  in  the  3-D  case). 


4  NUMERICAL  RESULTS  FOR  2-D  CASE 


The  two-dimensional  fourth  order  compact  scheme  is  applied  to  a  shock  reflection  problem  sketched 
in  Figure  1.  It  shows  a  5°  shock  at  Mach  1.95  reflecting  from  a  flat  plate.  Results  are  presented 
both  without  any  explicit  artiflcial  viscosity  (AV)  and  with  a  fourth  order  AV  term  of  the  type 
—  added  to  the  right  hand  side.  Addition  of  this  AV  term  reduced  the  accuracy  of  the 

competct  scheme  to  third  order  (note  that  in  3-D  the  sixth  order  AV  terms  precludes  this  reduction 
in  accuracy).  Calculations  are  also  made  with  a  second  order  implicit  Euler  AF  scheme  without  and 
with  the  preceding  AV  term.  Figures  2a  and  2b  show  the  results  of  the  compact  scheme  whereas 
Figures  3a  and  3b  show  the  corresponding  results  from  the  second  order  scheme.  It  is  seen  that  for 
e  =  0,  both  fourth  and  second  order  schemes  produce  very  oscillatory  results  but  with  e  =  0.36,  the 
restilts  from  the  comprct  scheme  (Figure  2b)  improve  dramatically  and,  in  fact,  are  much  better 
than  the  corresponding  results  obtained  from  the  second  order  scheme  (Figure  3b).  The  shocks 
captured  by  the  compact  scheme  are  sharper  and  the  convergence  is  also  seen  to  improve. 

Results  from  a  study  of  grid  aspect  ratio  effect  'tn  the  compact  scheme  are  also  presented  here 
in  terms  of  the  similarity  parameter  k  of  reference  1.  Three  values  of  k  are  considered,  namely  1.67, 
1.01,  and  0.42.  Figure  2b  corresponds  to  k  =  1.67  and  figures  4  and  5  correspond  to  k  =  1.01  and 
0.42,  respectively.  In  all  these  cases,  e  is  set  equal  to  0.36.  It  is  clear  from  the  figures  displaying 
the  effect  of  k  that  the  best  results  are  obtained  for  k  near  unity.  In  reference  1,  a  linear  theory 
(e.g.,  for  weak  shock)  predicts  the  same  results.  It  is  interesting  that  we  find  numerically  that  this 
is  also  the  case  in  the  present  nonlinear  problem. 

SXJMMARY 

1.  The  steady  state  solution  of  the  Euler  equations  of  gas  dynamics  may  be  achieved  to  fourth 
order  accuracy  using  a  compact  grid  stencil  of  3  x  3  and  3  x  3  x  3  in  the  2-D  and  3-D  cases 
respectively.  We  presented  two  examples  of  such  algorithms:  one  implicit  (Euler  approximate 
factorizations  scheme)  and  one  explicit  (Four-stage  Runge-Kutta). 

2.  Numerical  experiments  were  carried  out  for  the  2-D  shock  reflection  problem,  using  the  im¬ 
plicit  algorithm.  Comparisons  are  made  with  a  corresponding  second  order  scheme.  The 
results  show  that  the  compact  higher  order  scheme  offers  marked  improvement  in  both  accu¬ 
racy  and  convergence  rate. 

3.  In  connection  with  this  work  we  would  like  to  make  the  following  remarks.  It  is  known  that 
the  finite  difference  scheme  cannot  obtain  high  order  accuracy  for  conservation-form  equations 
computed  on  a  non-uniform  grid.  This  observation  coupled  with  the  ease  of  obtaining  fourth 
order  compact  schemes  even  in  3-D  for  the  Euler  equations  revives  an  old  debate:  Can  one  use 
uniform  grid  and  apply  conveniently  boundary  conditions  to  arbitrary  shapes.  The  potential 
gain  in  reduced  number  of  computational  nodes  and  enhanced  convergence  rate  appears  large 
enough  to  study  this  question  again. 
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Figure  1:  Shock  reflection  problem 


